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Abstract 

m 

Wc discuss steady viscous Kcplerian decretion discs around rapidly rotating stars. We assume that low 
^D I frequency modes, which may be excited by the opacity bump mechanism, convective motion in the core, or 

^^ ■ tidal force if the star is in a binary system, can transport an enough amount of angular momentum to the 

^ I region close to the stellar surface. Under this assumption, we construct a star-disc system, in which there 

forms a viscous decretion disc around a rapidly rotating star because of the angular momentum supply. 
We find a series of solutions of steady viscous decretion discs around a rapidly rotating star that extend 
^^ , to i?disc > WR* with i?* being the equatorial radius of the star, depending on the amount of angular 

O^ ' momentum supply. If the supply is low, the extension i?disc becomes small. 
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Introduction 
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1-G ] Discs around Be stars are now believed to be viscous Keplerian decretion discs (e.g.. Porter & Rivinius 2003; 

^~-*' Lee, Saio, & Osaki 1991). Formation mechanism for such decretion discs have not yet been identified, but it has been 
Q I suggested that angular momentum transferred by non-axisymmetric oscillation modes could play a role in the formation 
p) . of viscous discs around rapidly rotating B stars (e.g., Ando 1983, 1986; Lee & Saio 1993). Angular momentum transfer 
c/3 ' by non-axisymmetric oscillation modes in rotating stars has been discussed by various authors. Circularization of 
binary orbit and synchronization between the spin of a massive star and the orbital motion of the companion in a 
binary system are believed to result from angular momentum exchange between the orbital motion and the spin of 
the stars, where angular momentum redistribution in the star takes place through dissipative low frequency 5-modes 
tidally excited by the orbital motion of the companion (e.g., Zahn 1975, 1977; Goldreich &; Nicholson 1989; Papaloizou 
& Savonije 1997; Witte & Savonije 2001; Willems, van Hoolst, & Smeyers 2003). The problem of angular momentum 
transport in the Sun by low frequency g-modes has long been a hot topic in the field since the early years of helio- 
v^ ' seismology and investigated by many researchers including Press (1981), Schatzman (1993), Gough (1997), Kumar, 
Talon, & Zhan (1999), Talon, Kumar, & Zahn (2002), Mathis et al. (2008), where low frequency modes are assumed 
to be excited by convection motion in the convective envelope and suffer radiative damping as they propagate into 
the radiative core. 

To discuss angular momentum transfer by non-axisymmetric oscillations in a rotating star, we use a theory of wave- 
meanflow interaction, in which fluid motions are separated into waves and the mean flow and dissipative processes 
of the waves propagating in the mean flow have an essential role for the forcing on the meanflow (e.g., Andrews & 
Mclntyre 1978ab; Dunkerton 1980; Grimshaw 1984; Ando 1983; Goldreich & Nicholson 1989). For rotating stars, we 
Vh ■ regard the rotation velocity field as the meanflow and global oscillations as waves. In the Cowling approximation, in 
which the Euler perturbation of the gravitational potential is neglected, the meanflow equation may be given by (see 
the Appendix) 



p-4/ = yIm[V.(CV)], (1) 

where £ and p' are respectively the displacement vector and the Eulerian pressure perturbation of an oscillation mode 

and i = fsin^w^ is the specific angular momentum and (^) = (27i')^^ Jn idcfi denotes its zonal average, and we have 

separated the coordinates Xa into two parts such that Xa=Xa+ £,a (t,x), where Xa are the Lagrangian mean coordinates 
in the sense that (^q) = with ^a being the displacement vector. The right-hand-side of equation (1) represents the 
forcing effects associated with the oscillation mode, and an example of the forcing term calculated for low frequency 
modes in a slowly pulsating B star is given in the Appendix. 

In this paper, we try to construct a disc-star model, which is composed of a rapidly rotating star and a decretion 
disc around it. To support a decretion disc around a rotating star an enough amount of angular momentum must 
be supplied to the surface layers of the star, and this angular momentum supply is assumed to result from angular 
momentum deposition by axisymmetric oscillation modes. The rotation velocity in the surface layers is accelerated 
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acquiring angular momentum transferred by non-axisymmetric low frequency modes and a viscous Keplerian decretion 
disc forms around the star with the excessive angular momentum. Since the K-mechanism associated with the iron 
opacity bump excites low frequency and high radial order g-modcs and r-modes in slowly pulsating B (SPB) stars 
and low radial order /- and p-modes in (3 Ccphei stars (Dzicmbowski et al. 1993, Gautschy & Saio 1993), we may use 
the low frequency modes as an agent which brings about angular momentum transfer for SPBe stars. For early Be 
stars, however, the opacity bump mechanism does not work for driving low frequency modes, and we may resort to 
stochastic excitation of low frequency modes by convective motion in the core (Neiner et al 2012). If the Be star is in 
a binary system, low frequency modes may be tidally excited by the orbital motion of the companion star and they 
could work for disc formation. 

Okazaki (2001) has calculated steady and transonic decretion discs around Be stars, taking account of the effects of 
the radiation pressure on the radial flow. He has shown that the velocity field in the disc is very close to the Keplerian 
near the star and tends to angular momentum conserving in the region far from the star and that the effects of the 
radiation pressure are not essential for the existence of the transonic disc solutions. In his numerical analysis, he 
did not treat a rotating star and a disc around it as one system and assumed no angular momentum source in the 
subsurface layers of the star. These two points are what we concern in this paper, and we are interested in the part 
of the system around the rotating star. 

In the angular momentum conservation equation for a disc-star system, we need to consider both the forcing term, 
which appears in the meanflow equation and may be regarded as a source of angular momentum in the subsurface 
layers of the star, and the term due to the viscous torque, which is essential for angular momentum transfer in the disc. 
We derive in §2 a set of ordinary differential equations for a disc-star system, in which a Keplerian viscous decretion 
disc is assumed to form around a rapidly rotating star. In §3, we give numerical results obtained for steady disc-star 
systems, and we conclude in §4. In the Appendix we derive the meanflow equation we use in this paper. 

2. Equations for Decretion Discs around Rotating Stars 

For rapidly rotating SPB stars, for example, numerous prograde sectoral g-modes and r-modcs are destabilized 
by the opacity bump mechanism (e.g., Aprilia, Lee, & Saio 2012). Knowing their amplitudes and summing up all 
accelerating and decelerating contributions to the forcing on the meanflow, we can in principle estimate the amount 
of angular momentum deposited to the surface layers. With a linear theory of oscillations, however, we have no means 
to determine the amplitudes. In this paper, assuming that the contributions to accelerating the surface layers are 
dominant over those to deceleration, we consider a disc-star system, in which there forms around a rapidly rotating 
star a viscous Keplerian decretion disc, supported by excessive angular momentum transferred by the low frequency 
modes. 

To treat both a rotating star and a viscous disc around it as one system, we work in cylindrical coordinates {R,ip,z). 
The displacement vector of an oscillation mode may be given by ^ = CrGii + ^ipe^p +^z^z- If the amount of angular 
momentum (torque) provided by the low frequency modes is large enough to support a Keplerian viscous decretion 
disc around the star, the equation for angular momentum conservation, which includes both the viscous torque term 
for the disc and the forcing term due to the non-axisymmetric oscillations in the surface layers of the star, may be 
given by 



d 1 C^ / 2 N 771 



1 d , 9(Cp') 



R dR ' ^'^'^ ' dz 



(2) 



where aji^p is the Rip component of the viscous tensor and is assumed to dominate other components of the tensor for 
a geometrically thin disc. Here, we assume that the system is axisymmetric and in a steady state when the oscillation 
amplitudes are saturated, for example, by non-linear couplings between the oscillation modes (e.g., Lee 2012). If this 
is the case, integrating vertically equation (2), we obtain 

where we have used an^p = r]R{dfl/dR) with rj being the shear viscosity coefficient, v^ = RVl, and we have assumed 
p' = at the surface given by z = zq{R). Here, M = 2t: RY,Vfj is the mass decretion rate, which is assumed a constant, 
and E = J^° pdz is the surface density. Integrating equation (3) from _Ro to R, we obtain 

2^^'^/ vdz = giR)^MjiR)-fiR)~J, (4) 



dR.^z, 



where 

j(i?) = R^n, (5) 
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jmrR 



dz lm{CRp'), 



J = Mjn - ,90 - /o = Mj -g- f, 

and jo, go, and /o denote the quantities evaluated at R = Ro- 

Following Paezyrisky (1991), we also assume the disc-star system is in hydrostatic balance: 
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(6) 

(7) 

(8) 
(9) 



where the gravitational potential is given by $ = —GM^,| \/^F+z^ with M* being the mass of the star. Considering 
two neighboring points {R,zo) and {R + 6R,zo + 5zo) on the surface of the disc-star system such that p{R,zo) = = 
p{R + SR,zo + Szo), we obtain, using equations (8) and (9) for hydrostatic balance. 
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(10) 



dzo _ R fr^n"^ 
'dR~~o \GM^ 

where r = V-R^ + z'^. 

Rewriting equation (10) as dr^ /dR^ — r^fl/GM,, and assuming Q,/GMi, is a constant, we integrate equation (10) to 
obtain 



z^oiR) 



2GKh 



R' 



1/2 



(11) 



where Rp = Zo{Q) is the polar radius of the star. Using the condition Zo{Re) = 0, we may define, as a function of Rp and 
ri, the equatorial radius i?e of a star that is uniformly rotating at a rate Vt. We can also define the critical equatorial 
radius i?* of a star rotating at the critical rate ilo = (GM^/i?^)^/^. For this critical radius, we have R^/Rp = 1.5. 
Using this critical radius, we may rewrite (11) as 



yoix) 



1 



n^x^ 



1/2 



(12) 



where x = R/R^,, yo{x) = Zq{R)/R^,, yp = Rp/R^, = 2/3, and J7 = fl/rto- Note that yo{x) = gives the critical radius 
x = 1 at rj = 1 for a uniformly rotating star. The rotation rate fl may be normalized by J7e = (GM^/i^g)^/^, and in 

this case we have fi/fie = Xe fi with Xe = Re/ R*, where Xe is determined as a solution to yo{x) ~ for a given fj < 1. 
Figure 1 plots Xe and O/rie as a function of f2 = Vt/Vlo- As fj decreases, Xe tends to 2/3. We also note the rapid 
decrease of ^/VL^ as O decreases from 1. 

We divide the disc-star system into two parts, that is, the inner part (i? < i?tr) and the outer part {R > i?tr)- We 
assume the inner part is uniformly rotating at a constant rate Jig and the surface of the inner part is given by equation 
(11) with D, = Vis- On the other hand, the outer part of the system is composed of the outer part of the rotating star 
and a decretion disc and is allowed to rotate differentially. We assume the outer part extends to the radius i?out 3> R* ■ 

For a thin disc, the shear viscosity coefficient rj may be given as 

'q = apoVsfiZo, (13) 

where a is a dimensionless constant parameter such that < a < 1, and po and Vs,o = {dPo/dpoY^^ are the density 
and the sound speed evaluated at the equatorial plane. Assuming a polytropic relation P = Kp^~^^'" , wc integrate 
equation (9) to obtain 



P = Pa 



and hence 
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Fig. 1. Xe (solid line) and C/Qe (dashed line) as a function of f^ = C/Co, whore Qe = {GhU/Rlf/'^ and f^o = {GM,/Rl)^^^, and 
Re and _R» arc the equatorial radius of a star uniformly rotating at Q and Qq, respectively, and M, denotes the mass of the star. 
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Equation (4) is now given by 



9' 



Ana 



{Ghh 



in+0.5 



XT. 



^2n+3 
^0 



da 



i?3"-l-5 dR 



MR^n-f-J. 



(16) 



(17) 



Equations (10) and (17) make a set of ordinary differential equations we have to solve with appropriate boundary 
conditions imposed at i? = i?tr and R = Rout ■ 

In this paper, we assume n= 1.5 for the polytropic index. Equations (10) and (17) then reduce to 

dy 



X 

dx y 

dn _lx^ - 
dx a 2/* 
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(18) 
(19) 

(20) 

(21) 



and p = A/,/(47ri?2/3), and the summation in equation (21) is over the oscillation modes which contribute to the 
forcing on the meanflow, where k denotes a collective mode index. For the polytropic index n — 1.5, we may use 
K ^ 0.4242GA'/y^i?* (Chandrasekhar 1939) to obtain 

1/2 

(22) 



M, 



a ~ Ana ( — r^ I [ * 

and 

po ~ 1.5p y^x''/^ 



(23) 
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Although the constant a depends on the viscosity parameter a, the magnitude of a is almost the same as that of the 
constant b, that is, a ^ b. As suggested by Paczyhsky (1991), the constant a can be as large as a ^ lO^'^a, depending 
on the quantities such as M, M^, and i?*. As a increases, however, it becomes difficult to numerically find solutions 
to the set of differential equations. In this paper, we employ a = 10^, which leads to the disc thickness zq/R^ 0.1 (see 
below). 

To determine the function f{x) with numerous low frequency non-adiabatic modes of a rapidly rotating star, we 
need to correctly estimate their amplitudes, which are difficult with a linear theory of oscillation. For our numerical 
computation, assuming that the resulting forcing function f{x) is that for accelerating the rotation velocity in the 
surface layers, we use a simple function given by 

bf{x)=^bo( p-— ^ ^— -l)^5o/o(x), (24) 

\exp[-bi[x- xoXe)\ + l ) 

where 6o: ^ij and xq are parameters, and x^ ~ Re/ R*- In the surface region of the star where a; ~ Xe < 1, we may 

approximate |/(xi = a;oa;e)|-a;^z/o(a;i)(/Oo/p)|Efc'^Im[(C^A)(p'/A90] I ^a^r^'^yo(2;i)IEfe'^/|sin4| where A- - |^fl,/r| 
is the normalized oscillation amplitude and 5k is the phase difference between (^n and p' near the surface, and we have 
approximated gr/{GM<,/R<,) ~ 1. Since |/o(a;i)| = 0.5, we may have &o ~ ^ x Xi'^'^yoi^i)\J2k^^fk^™^k\- Although it 
is difficult to correctly estimate the magnitudes of the quantities such as |X)fc"^/fc siii'^fcl, assuming |Ej,msinJfe| ^ 1 
and Xi^'^yQ{xi) ~ 10~^, we have 6o ^ 1 for fk ^ 10"^ and b^ a = 10^. Note that for a: ^ 1, we should have bf{x) = 0, 
expecting there occurs no forcing on the meanflow in the disc. 

The set of differential equations (18) and (19) are integrated with three boundary conditions, one given by dCl/dx — 
dClx/dx with i),K = x~^''^ at a; = Xout = Rout/R*, and other two conditions by ^{xtr) = ^s and y(xtr) = yo(a;tr) 
at x = xtr = Rti/R*- Three boundary conditions for the second order ordinary differential equation are needed to 
determine the parameter j. In this paper, we use xtr — 0-8. 

As indicated by equation (19), for x 3> 1, dCl/dx changes its sign at Xj, for which x^Cl(xj) — J = 0, where we have 

assumed bf{x) — 0. If we assume ft ex CIk, we obtain Xj ~ j^. To understand a rough property of the solution in the 
region l<$i x <$iXj in which x^Cl'^ j, assuming functional forms given by y =px^ and J7 = qCIk, where the parameters 
p, q, and s are assumed only weakly dependent on x, and substituting the expressions into equations (18) and (19), 
we obtain 

11 / ■■ \ ^^^ 1 

^ = T^' ^H4j ' '^^(l+,..-Va)V- (25) 

These relations may be consistently satisfied \i p ^ 0.1 and q ^ 1 for j ~ 10 and a ~ 10^, which indicates that the 
decretion discs are geometrically thin for the parameter values. If the disc extends beyond x^, it may become angular 
momentum conserving so that x^Ct — j be a constant (see Okazaki 2001). If this is the case, the outer boundary 
condition must be modified and an appropriate treatment of solutions around a point of Xj will be required. 

3. Numerical Results 

For the case of Os > 1, we integrate equations (18) and (19) using a Henyey type relaxation method with the initial 
guesses given by y = a;[(jx~^/^ — l)/1.5a]^/^ and fl = x~^^'^, and find the value of j that satisfies the three boundary 
conditions. For the case of Clg < 1, since there appears near a; ~ 1 a point Xm, at which dVt/dx = 0, we integrate, for 
a given value of j, the differential equations (18) and (19), with an implicit Rungc-Kutta method, from x = a^tr using 
y{xf[) = yo{xf[) and Vl{xfc) = ^s, to tlie point Xm- Since this integration gives Xm, ym = y{xm), and Ctm = ^{xm) as 
a function of j , or equivalently, y^n, f^m, and j as a function of Xm, we solve equations (18) and (19) in the interval 
between Xm and Xout with the relaxation method using the initial guesses to find the value of Xm such that the 
boundary conditions y(xm) = y-m and Q,{xm) = ^m at x = Xm and dCl/dx ~ dClx/dx at x = Xout are satisfied. 

Figure 2 shows y and Jl as a function of x for Cls = 0.98, 0.99, 1.00, and 1.01, where we have assumed 6o = 6i = 50, 
a;out = 10, and xq = 0.95. In Table 1, we tabulate several characteristic quantities such as Xm and j as a function of 
rJs. As shown by the left panel of the figure, there appears a sharp dip in y at the boundary between the star and disc 
when n^ < 1, and the dip becomes deeper for smaller values of Cls- Tabic 1 indicates that this star-disc boundary is 
located at a radius near Xe for il^ < 1. For the case of fJs <1, if we go outwards from xtr, at a point near Xe fl starts to 
rapidly increase to attain a super-Keplerian rate {ft > Hk) at Xm and then decreases to the Keplerian velocity. Note 
that f2 is slightly sub-Keplerian in the region of x >1. However, for 51^ = 1.01, we have Q < Qk in the entire x region 

of the solution and there appear no points of dQ/dx = 0. For a given bg, there exists the lower limit of ilg, below 
which no solutions to the differential equations are found. As Og decreases from 1 , the amount of angular momentum 
deposition required to accelerate the sub-Keplcrian rotation velocity to a super-Kcplerian one is increased, and hence 
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Table 1. Xm ■■ 


and j for dccrotion disc solutions 


ns/iio 


Xe 


iis/rte 


'^rn 


3 




bo = 50 


xo = 0.95 


.Tout = 10 




0.9750 


0.890 


0.819 


0.893 


5.211 


0.9800 


0.900 


0.837 


0.905 


4.859 


0.9850 


0.912 


0.858 


0.919 


4.405 


0.9900 


0.926 


0.883 


0.936 


3.837 


0.9950 


0.946 


0.916 


0.957 


3.603 


0.9980 


0.965 


0.946 


0.963 


5.420 


0.9990 


0.975 


0.962 


0.964 


7.629 


1.0000 


1.000 


1.000 


0.965 


16.81 


1.0100 








52.98 




6o = 10 


a;o = 0.95 


3^out — -^ 




0.9935 


0.940 


0.905 


0.945 


1.688 


0.9940 


0.942 


0.908 


0.948 


1.666 


0.9950 


0.946 


0.916 


0.953 


1.619 


0.9960 


0.952 


0.925 


0.960 


1.578 


0.9970 


0.958 


0.935 


0.966 


1.552 


0.9980 


0.965 


0.946 


0.972 


1.577 


0.9990 


0.975 


0.962 


0.976 


1.761 


1.0000 


1.000 


1.000 


0.976 


3.160 


1.0100 








37.33 




&o = 50 


a:o = 0.99 


a^out = 10 




0.9500 


0.854 


0.750 


0.856 


19.41 


0.9600 


0.867 


0.774 


0.870 


18.77 


0.9700 


0.881 


0.803 


0.886 


17.82 


0.9800 


0.900 


0.837 


0.907 


16.39 


0.9900 


0.926 


0.883 


0.936 


15.10 


0.9950 


0.946 


0.916 


0.952 


16.86 


0.9980 


0.965 


0.946 


0.962 


21.97 


0.9990 


0.975 


0.962 


0.966 


25.70 


1.0000 


1.000 


1.000 


0.973 


36.16 


1.0100 








66.14 



the derivative d£L/dx inevitably becomes steeper in the region where acceleration takes place. The lower limit of £7^ 
is determined by the value at which dCl/dx = oo or y = 0. For 6o = 50 and Xout = 10, the lower limit of Cls is ~ 0.975. 

Figure 3 shows that the ratio y/x is less than ~ 0.1 in the disc, indicating the disc is geometrically thin. Since the 
ratio is approximately proportional to j^/^ as suggested by equation (25) and the value of j for Vlg = 0.98 is larger 
than that for Vis = 0.99 (see Tabic 1), the ratio y/x for the former is larger than that for the latter. This figure also 
shows that the ratio gradually decreases as x increases from x ~ 1 . 

As suggested by Table 1, proper solutions to the differential equations can be obtained only when the outer boundary 
condition is imposed at a;out;< Xj ~ j^, and the properties of the solutions for a given j do not strongly depend on 
Xout so long as Xout<Xj for Xout > Xm- If we employ &o = 10 instead of &o = 50, we can find solutions for Xout = 2 but 
no solutions for .Tout = 10 since the parameter value of j we obtain for 6o = 10 is ^ 1.5 for il^ < 1. The lower limit of 
fJs for 6o = 10 is ~ 0.9935, which is much closer to unity than the lower limit 0.975 for the case of h^ = 50. Figure 4 
plots y and H as a function of x for 6o = 10. The discs for b^ = 10 are thinner than for &o = 50, and the peak value of 
VL approximately corresponding to the lower limit of fJs is smaller than that for the case of &o = 50. 

As ^s approaches unity from $7^ < 1, the parameter value j begins to increase rapidly (see Table 1). The large j 
means that decretion discs having large radii are possible, for example, a disc can have an extension of Tout >100 if 

j > 10. 

Assuming that the acceleration takes place in the region very close to the stellar surface, that is, in the region of 
much lower density, we examine the case of Tq = 0.99, where we use 6o = 6i = 50 and Tout = 10. Wc obtain a series of 
decretion disc solutions, the properties of which are quite similar to those for tq = 0.95, except for that the values of 
j for To = 0.99 are greater than those for tq = 0.95. Because of the large values of the parameter j, the discs can have 
a larger extension for to = 0.99 than for tq = 0.95. The lower limit of Clg is ~ 0.95, which is smaller than ~ 0.975 for 
the case of tq = 0.95. 
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y 0.4 
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Fig. 2. y = zo(R)/R, (left panel) and Q, = Q/Qq (right panel) as a function of a; = R/R, for Qa = 0.98, 0.99, 1.00, and 1.01 for 
hQ = b\ = 50, xq = 0.95, and Xout = 10, where f2o = (GMt/R^)^'^ with M» being the mass of the star, and fls denotes the value of 
n at X = xtr = 0.8. The dotted line in the right panel indicates the Keplerian rotation rate given by C^ = x^'"^'^. 
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Fig. 3. The ratio y/x versus x = R/R^ for Q^ = 0.98, 0.99, 1.00, and 1.01 for bo = bi = 50 and xq = 0.95. 
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Fig. 4. y = zo{R)/R, (left panel) and Q, = U/flo (right panel) as a function of x = R/R^, for fls = 0.994, 0.996, 0.998, 1.00, and 
1.01 for bo = 10, bi = 50, xq = 0.95, and Xout = 2. The dotted line in the right panel indicates the Keplerian rotation rate given by 
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Fig. 5. y = zo{R)/R, (left panel) and C = U/Uo (right panel) as a function of x = R/R, for Us = 0.95, 0.96, 0.97, 0.98, 0.99, 1.00, 
and 1.01 for feo = ''i = 50, Xout = 10. and xq = 0.99. The dotted lino in the right panel indicates the Keplerian rotation rate given by 
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4. Conclusion 

Wc have obtained a series of solutions for steady viscous Keplerian decrction discs around a rapidly rotating star, 
assuming the existence of angular momentum supply to the region close to the surface of the star. The angular momen- 
tum supply may be provided by angular momentum deposition that takes place through wave-meanflow interaction, 
where the waves are low frequency global oscillations excited by the opacity bump mechanism for SPB star, or by a 
stochastic mechanism for early Be stars, or by a tidal force if the star is in a binary system. We may conclude that 
the angular momentum supply to the surface layers can be a mechanism for disc formation around rapidly rotating 
Be stars. 

If the amount of angular momentum supply, which may be represented by the parameter bo in this paper, is large 
enough, viscous decrction discs can extend to a distance as far as i?out > lOi?*, but the possible disc radius may 
decrease as the supply becomes small. We also find that if the acceleration occurs in the region very close to the 
stellar surface, the possible extension a decrction disc attain can be as large as i?out >100-R». For a given bo, there 
exists the lower limit of Og, below which no steady decrction disc solutions are found. If bo is large enough, decrction 
disc solutions can be obtained even for the rotation rate ^--^80% of the critical rate fie determined by the equatorial 
radius Re of the star, and the disc radius can be large enough for a Be star. If bo is small, however, disc solutions can 
be found only when Cls is very close to unity. If the angular momentum supply is provided by global oscillations, the 
parameter Bq is closely related to the amplitudes of the oscillations. We have argued that the amplitudes of order of 
^r/R* ^ 0.01 — 0.1 can lead to reasonable values of &o- 

The density at the mid-plane of the disc may be estimated by using equation (23), which leads to po oc x~^'^ if the 
ratio y/x is assumed almost constant. Since po ex x~^'^ has been suggested observationally, the x dependence of po in 
our model is in a serious conflict with the observational estimation (e.g.. Porter & Rivinius 2003). For B type stars 
we have p ~ 10~^, and since y/x ~ 0.01 — 0.1, we obtain po ~ 10~^ — 10~^ at x ~ 10, the value of which is much higher 
than that observationally estimated (e.g.. Waters 1986). We may be able to use much larger values for the parameter 
a to reduce the ratio y/x, but it becomes numerically difficult to obtain disc solutions to the differential equations. 

As indicated by the plots of ^{R), there occurs a strong differential rotation in the region close to the surface, 
particularly for lower values of fJs . The strong differential rotation can change the modal properties of oscillations and 
hence the contributing to the forcing on the velocity field. Stability analysis of low frequency modes in differentially 
rotating star, which will be one of our future studies, is necessary if for the forcing we resort to the modes that are 
excited by the opacity bump mechanism. 

Appendix 1. Mean Flow Equation in the Lagrangian Mean Formalism 

Following Grimshaw (1984), in this Appendix we derive a meanflow equation for zonal flows around the rotation 
axis of stars, using the Lagrangian mean formalism (see also Andrews & Mclntyre 1978b). In a frame rotating with 
the angular velocity Jlc, the cj) component of the momentum conservation equation in spherical polar coordinates 

x={xi,X2,x:>,)=[r,9,(j)\ is given by 



dv^, VrV^ , weu^cotfl , „o / • a- , a- 
at r r V 



I dp „ 1 a* ,,^, 

/'T-^— -' (Al) 



' d<j) fsmO d(j) 



df ^ ^dO ^ ^ . ;:d4> 

ve = r—, v^^rsmO—, (A2) 



where 

and 

dt dt dr r QO fsmedcj) 

To discuss wave- mean flow interactions, wc introduce the Lagrangian mean coordinates x = {xi,X2,X3) = {r,9,(j)) and 
the displacement vector ^a for a = 1, 2, 3 such that 

Xa^Xa+£,a{t,x). (A4) 

We assume that for any given Va there is a unique "reference" velocity v^, such that when the point x^ moves with 
velocity Va the point Xq, moves with velocity Va (e.g., Grimshaw 1984). Using the reference velocity, we may define 

1 = 1+^,^A + ^A + ^^A, (A5) 

dt dt ^ dr r 36 rsinO dd>' 
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where d/dt indicates the partial time derivative with the coordinates x being held constant. For the wave-mean flow 
interaction formulatoin for zonal flows, we introduce an averaging procedure defined by 

{f)^^rfd^, (A6) 

and wc assume for the displacement ^q, 

(ea)=0, (A7) 

where (p = x^. The velocity v may be regarded as the mean velocity associated with the coordinates {xa)^ and the 
displacement ^q, represents the waves. 

Using the Jacobian J for the coordinate transformation between {xa) and {xa) given by 

we define the mean density p associated with the coordinate x as 

pf^sin^J = pr^sin6'. (A9) 

It is convenient to introduce (e.g., Andrews & Mclntyrc 1978b) 

for which we have 

^=0, (All) 

dxp 



and 



\J iXj (y^ ^ yj tLi ly 



i^«^^=^/3.J = i^/3a^, (A12) 

where we have employed the dual summation convention that Greek indices are summed over the range 1 to 3, and 
£123 = £231 = £312 = 1 = -£i32 = -£213 = "£321 and £0/37 = Otherwise. 
We rewrite equation (Al) as 

.dl d ^ ^d ^ ,, ^, 

where £ is the specific angular momentum, in an inertial frame, around the rotation axis defined by 

i = fam§v^ + flcr^sm^§. (A14) 

Multiplying equation (A13) by f^sin^J, we obtain 

di dp 1 d d ^ ,. , 

dt ocp r'^smOoxp d4> 



s2„:„fl^ „2-„a^\ ^2„-n'-^^l 



where 

-R3/3 = (^3/3 ( Jr^ sin ^p-r^ sin 6*^1 - f'^ sinOp^^K^iJ, (A16) 

and p ~ p(p,s) with s being the specific entropy, and we have used the identity 

K^pr^ smOp = d^liJf'^ sinOp ~ r^ sinep-^K^p. (A17) 

Applying the averaging procedure (A6) to equation (A15), we get 
dll) I d .„ ^ /^d 



dt r'^awwoxp \ dcf) / 

where we have used {dp/ dcf)) ~ 0. In general, 

;;2,.- fl- ^2„- /)--\ / ^'2 „■ a^'^^l 



(Rap) = S„p{jf'^sm0p-r^smep) - ( f^sm0p-^K^0 ) (A19) 
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is called the radiation stress tensor (e.g., Grimshaw 1984). 

Equation (A18) may be regarded as the (f) component of the meanflow equation, the left-hand-side of which may 
represent the time evolution of the meanflow and the right-hand-side the forcing by the waves represented by S,a ■ So far 
we have not assumed that the amplitudes of the waves S,a are infinitesimally small, and in principle we can formulate 
the wave-meanflow interaction as a nonlinear theory, which includes equations of motion for both the meanflows and 
waves. To avoid solving such a difficult non-linear problem, we use a linear theory to describe waves S,a a-nd we 
calculate the forcing terms in the meanflow equation using the linear waves ^a ■ 

If we employ a linear theory to describe waves represented by the displacement ^q, , we may write 

Of) 
p^p + 5p = p + p' + U^, (A20) 

where p is the pressure in equilibrium state, 5p and p' denote the Lagrangian and Eulcrian perturbations, respectively. 
Applying the averaging procedure (A6), to second order of the perturbations we obtain after some manipulations for 
the second term on the right hand side of equation (A16) 



'''''''•>iSH'\''"'"'sS')^\s:^[''''"9!^J^'' >■ '*'" 



where we have used [df/dcj)) ~ d (/) / dcj) = and dp/d<j) = for the equilibrium pressure p, and x^ = 4>. Because 
pr smt'-— — 4^, 



^^^fpr^sine^A^ — 
dxa dxy \ dcj) ''I d(f> 






\J tij Q/^ kJJj'^ 
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pr^sinei^^^ ), (A22) 



d d ( ^ ^ d^^ 



we 



find 



d 5 ^^2,. A,\\_/ d a /_2 . „, d^^ 



^pVsin.^^.J ; - \a^a^ [P^'^-^-'^-H) ) - 0' (^23) 

and hence we can omit the second term on the right-hand-side of equation (A21) to obtain 

'f^ sinep^K^o) = (r^ sin0^p'\ , (A24) 



^2„-, fl-: „2„-, /)--\ / „2„- /jC'Ca^/ 



and hence 

{R3a) = ha Ur^ sin9p - r^ sui9p) - (r^ sind^p' ) . (A25) 

The meanflow equation (A18) is then reduced to 



where ^ = ^r-Er + ^9^8 + ^0^0 with e,., eg, and e^ being the orthonormal vectors in the r, 0, and (p directions, and 
S,r = ^i;£,e = r£,2, and ^0 = r sin 6*^3, and p = p. Since 

which is correct to second order of perturbations, we obtain 




where we have used d^/dcj) ~ 0, V^<f>' = inGp' and (V$' • (9V<f>'/90) = 0. Since the azimuthal and temporal depen- 
dence of the perturbations are assumed to be given by the factor exp (Ittk/) -|- iujt) with m and to being the azimuthal 
wavenumber and oscillation frequency, using, for example, 

iU^r) = ^Re (qe.) = l^ei^^C), (A29) 

we can rewrite equation (A28) as 



12 



d(i 



dt 



-Re 



-Im 



.* dp' , a$' v$'* a$' 

^ 'd^'^^^ 'd^ ^ 47rG ^ 
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(A30) 



where the asterisk indicates the complex conjugation. 
Integrating over a spherical surface, we obtain 



d{n m d , , 



pC*' 



ar 47rG 



(A31) 



This expression (A31) is essentially the same as that used by Papaloizou & Savonije (1997) who discussed the forcing 
by low frequency modes tidally excited in a massive star by the orbital motion of the companion in a binary system 
(see also Ryu & Goodman 1992, Lin, Papaloizou, & Kley 1993). In the Cowling approximation, we obtain 



dW{r) 



dt 
where 

W{r) 



2'Kr'^ dr 



(A32) 



Trr^Im (^7^ (A33) 

may be regarded as a work function (e.g., Unno et al. 1989). We note that dW/dr > and dW/dr < respectively 
indicate the excitation and damping regions for the oscillation modes. For uniformly rotating stars, non-axisymmetric 
(to 7^ 0) oscillations of rotating stars are separated into prograde and retrograde modes, and in our convention, positive 
(negative) m is used for retrograde (prograde) modes. From equation (A34), we find that there occurs acceleration 
(deceleration) of the zonal flow in the damping (excitation) regions of prograde modes {m < 0), while deceleration 
(acceleration) occurs in the damping (excitation) regions of retrograde modes (m > 0). 

It is possible to rewrite the mean flow equation (A26) using the Reynolds stress. Substituting into equation (A26) 
the component of the linearized momentum equation given by 



dp' 



-prsmt 



d_ 
di 



d 



where v^ — rsin0J7i(r, ^ 



tlO dcj) 



1 di 
rsinO dr 



1 ide 



rsin^ r dO 
and using 



ve- 



1 a$' 

iVL\0 d<f> 



7;' = -i+Oi^^e,-rsin0(^.Vr!i)e^ 



dt 



(A34) 



(A35) 



we obtain after some manipulations 



d(£ 



dt 



J_d_ 

j,2 Qj, 



{pr^sm9v',{v'^ + 2nfecoseCe)) 

Id I (9$' \ 3 

{prsin^Ov'g {v'^ + 2nframd(r)) - (p'^ > + 7^^ ' (P^sin^g^) 



rsin^ do 



dt 



dt \r'^ dr^ 



'■en Mr) 



1 



d 



rsmO do 



(prsin^ cos on fe£.g) 



where 

fr = 



(A36) 

(A37) 

(A38) 

If we neglect the terms like d{Q^) /dt and d(^^^) /dt assuming the amplitudes of the perturbations are saturated, for 
example, by nonlinear effects, the expression for the forcing terms in the mean flow equation (A36) with /e = 1 reduce 
to that proposed by Pantillon et al. (2007) and Mathis (2009). 

It is instructive to give an example of the forcing term in the meanflow equation. To see the behavior of the forcing 
term, we introduce the local timescale t^^ defined by 



lc)hir2 _ 1 dhin 

2 91nr ' '''' ^ "^ ^ 2 91nsin6'' 
I- 2nfe cos 0^0 + 2nfr sin O^r ■ 



pAM 



d{£)/dt 



(A39) 
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Fig. 6. 1/t for unstable I = \m\ = 1 322 modes for fic/{GM/ B'^)^'^ = 0.1 (left panel) and an unstable r^o mode for 

Uc/iGM/R^y^^ = 0.4 (right panel) of a AMq main sequence model with X = 0.7, Z = 0.02, and Xc = 0.3646, where the solid 
and dotted lines indicate retrograde and prograde modes, respectively. Note that uniform rotation is assumed. The amplitude 
normalization is given by Si^ (R) = 1 at the surface, and t^^ is given in units of seconds. 



An example of l/r is given for low frequency modes of a 4A-Iq main sequence star with X = 0.7 and Z = 0.02 in 
Figure 6, where the model has been calculated with a standard stellar evolution code with the OPAL opacity (Iglesias 
& Rogers 1996), and we have used the method of calculation given by Lee & Saio (1993) for non-adiabatic oscillation 
modes of a uniformly rotating star. As the figure shows, there occurs a strong acceleration by retrograde g- and 
r-modes in the layers at r/R ~ 0.95, although the prograde g-mode contributes to deceleration of the rotation. 
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